#include "TriMesh.h"
#include <map>
#include <queue>
#include "GeometryTools.h"

TriMesh::TriMesh(void)
{
}

TriMesh::~TriMesh(void)
{
}

void TriMesh::needBoundingBox()
{
	TriMesh::ConstVertexIter vIt = vertices_begin();
	TriMesh::ConstVertexIter vEnd = vertices_end();

    bbox_min = bbox_max = point(vIt);

    for (; vIt != vEnd; ++vIt) {
		const Point& p = point(vIt);
		bbox_min.minimize(p);
		bbox_max.maximize(p);
    }
}

void TriMesh::request_curvature()
{
	OpenMesh::VPropHandleT<TriMesh::Point> curvature;
	get_property_handle(curvature, "curvature");
	TriMesh::VertexIter v_it, v_end(vertices_end());
	TriMesh::VertexVertexIter vv_it, vv_it2;
	TriMesh::Scalar voronoiArea;
	for (v_it=vertices_begin(); v_it!=v_end; ++v_it)
	{      
		property(curvature,v_it)[0] = 0;
		property(curvature,v_it)[1] = 0;
		property(curvature,v_it)[2] = 0;
		voronoiArea = 0;
		for (vv_it=vv_iter( v_it ); vv_it; ++vv_it)
		{
			//         (p)  -  p3
			//        /    \     /  
			//      p1 -  p2
			// p is v_it; p2 is vv_it;
			Point p = point(v_it);
			Point p2 = point(vv_it);
			vv_it--;
			Point p3 = point(vv_it);
			vv_it++;
			vv_it++;
			Point p1 = point(vv_it);
			vv_it--;
			TriMesh::Scalar wij = ((p-p1)|(p2-p1)) / ((p-p1)%(p2-p1)).length() + 
				((p-p3)|(p2-p3)) / ((p-p3)%(p2-p3)).length();
			property(curvature,v_it) += wij * (p2 - p);
			voronoiArea += wij * (p2-p).length() * (p2-p).length() / 8;
		}
		property(curvature,v_it) /= 2*voronoiArea;
	}
}

void TriMesh::request_curvature_color()
{
	OpenMesh::VPropHandleT<TriMesh::Point> curvature;
	get_property_handle(curvature, "curvature");
	OpenMesh::VPropHandleT<TriMesh::Point> curvature_color;
	get_property_handle(curvature_color, "curvature_color");
	TriMesh::VertexIter v_it, v_end(vertices_end());
	std::vector<double> lens;
	TriMesh::Scalar len;
	for(v_it=vertices_begin(); v_it!=v_end; ++v_it)
	{
		len = property(curvature, v_it).length();
		lens.push_back(len);
	}
	std::sort(lens.begin(), lens.end());
	TriMesh::Scalar threshold = lens[ceil(9*n_vertices()/10.0)];
	for(v_it=vertices_begin(); v_it!=v_end; ++v_it)
	{
		if(property(curvature, v_it).length() > threshold)
			len = 1;
		else
			len = 0;
		property(curvature_color, v_it)[0] = 1-len;
		property(curvature_color, v_it)[1] = 1-len;
		property(curvature_color, v_it)[2] = 1-len;
	}
	lens.clear();
}

TriMesh* TriMesh::read(const char* filename, OpenMesh::IO::Options* opt)
{
	TriMesh* mesh = new TriMesh;

	OpenMesh::IO::Options default_opt;
	if (!opt)
		opt = &default_opt;

	if (!OpenMesh::IO::read_mesh(*mesh, filename, *opt)) {
		delete mesh;
		return NULL;
	}

	// If the file did not provide vertex normals, then calculate them
	if ( !opt->check(OpenMesh::IO::Options::VertexNormal) &&
		mesh->has_face_normals() && mesh->has_vertex_normals() ) {
		// let the mesh update the normals
		mesh->update_normals();
	}
	// add additional properties
	OpenMesh::VPropHandleT<TriMesh::Point> curvature;
	mesh->add_property(curvature, "curvature");
	OpenMesh::VPropHandleT<TriMesh::Point> curvature_color;
	mesh->add_property(curvature_color, "curvature_color");

	return mesh;
}

bool TriMesh::save(const char* filename,  TriMesh* mesh, OpenMesh::IO::Options* opt)
{
	OpenMesh::IO::Options default_opt;
	if (!opt)
		opt = &default_opt;

	if (!OpenMesh::IO::write_mesh(*mesh, filename, *opt)) 
		return false;
	return true;
}

void TriMesh::getMeshStat(std::vector<double> &meshStat, TriMesh* meshCopy)
{
	// statistics 1: percentage of vertices having 6 neighbors. [0]
	meshStat.clear();
	TriMesh::VertexIter v_it;
	TriMesh::VertexVertexIter vv_it;
	int valence;
	int n_6n = 0; // number of vertices having 6 neighbors
	for(v_it=vertices_begin(); v_it!=vertices_end(); ++v_it)
	{
		valence = 0;
		for(vv_it=vv_iter(v_it); vv_it; ++vv_it)
			valence++;
		if(valence == 6)
			n_6n++;
	}
	meshStat.push_back(1.0*n_6n/n_vertices());
	// statistics 2: radius ratio, 2*r/R. [1] Tmin, [2] Tmean
	TriMesh::FaceIter f_it;
	std::vector<double> radiusRatio;
	for(f_it=faces_begin(); f_it!=faces_end(); ++f_it)
	{
		radiusRatio.push_back(2 * GeometryTools::radiusOfInscribedCircle(this, f_it.handle()) 
			/ GeometryTools::radiusOfCircumscribedCircle(this, f_it.handle()));
	}
	meshStat.push_back(*std::min_element(radiusRatio.begin(), radiusRatio.end()));
	meshStat.push_back(std::accumulate(radiusRatio.begin(), radiusRatio.end(), 0.0)/n_faces());
	radiusRatio.clear();
	// statistics 3: min angles in triangles, [3] minimum degree, [4] mean degree
	std::vector<double> angle;
	for(f_it=faces_begin(); f_it!=faces_end(); ++f_it)
		angle.push_back(GeometryTools::minAngleInTriangle(this, f_it.handle()));
	meshStat.push_back(*std::min_element(angle.begin(), angle.end()) * 180.0 / 3.1415926);
	meshStat.push_back(std::accumulate(angle.begin(), angle.end(), 0.0)/n_faces() * 180.0 / 3.1415926);
	angle.clear();
	// statistics 4: relative RMS of edge lengths, [5] edge RRMS
	TriMesh::HalfedgeIter he_it, he_it_copy;
	double edgeRMS = 0;
	int num = 0;
	for(he_it=halfedges_begin(), he_it_copy=meshCopy->halfedges_begin(); he_it!=halfedges_end(); ++he_it, ++he_it_copy)
	{
		num++;
		TriMesh::Point p1_copy = meshCopy->point(meshCopy->from_vertex_handle(he_it_copy.handle()));
		TriMesh::Point p2_copy = meshCopy->point(meshCopy->to_vertex_handle(he_it_copy.handle()));
		TriMesh::Point p1 = point(from_vertex_handle(he_it.handle()));
		TriMesh::Point p2 = point(to_vertex_handle(he_it.handle()));
		edgeRMS += pow(abs((p1_copy-p2_copy).length() - (p1-p2).length())/(p1_copy-p2_copy).length(), 2);
	}
	edgeRMS = sqrt(edgeRMS / num);
	meshStat.push_back(edgeRMS);
	// statistics 5: relative RMS of volume magnitudes, [6]
	double volumeRMS = 0;
	double volume = 0;
	double volumeCopy = 0;
	TriMesh::FaceIter f_it_copy;
	TriMesh::FaceVertexIter fv_it, fv_it_copy;
	for(f_it=faces_begin(), f_it_copy=meshCopy->faces_begin(); f_it!=faces_end(); ++f_it, ++f_it_copy)
	{
		fv_it=fv_iter(f_it);
		TriMesh::Point p1 = point(fv_it.handle());
		++fv_it;
		TriMesh::Point p2 = point(fv_it.handle());
		++fv_it;
		TriMesh::Point p3 = point(fv_it.handle());
		volume += p3 | (p1%p2);
		fv_it_copy=meshCopy->fv_iter(f_it_copy);
		p1 = meshCopy->point(fv_it_copy.handle());
		++fv_it_copy;
		p2 = meshCopy->point(fv_it_copy.handle());
		++fv_it_copy;
		p3 = meshCopy->point(fv_it_copy.handle());
		volumeCopy += p3 | (p1%p2);
	}
	volumeRMS = abs(volume - volumeCopy)/volume;
	meshStat.push_back(volumeRMS);
}
